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The mult i- algorithm collaborative crystal structure prediction was realized and constructed as 
a Multi-algorithm-collaborative Universal Structure-prediction Environment (Muse). The evolu- 
tionary algorithm was coupled with the simulated annealing and the basin hopping algorithms to 
efficiently find the stable and met ast able structures of materials under certain conditions. After in- 
troduced two new operators, slip and twist, the diversity of structures are exceedingly enhanced. In 
particular, in order to achieve the self-adaptive evolution of structures, Muse uses the competition 
techniques of ten variation operators to increase the diversity of crystal structures. The multi- 
algorithm coupling, the ten variation operators, the symmetry constraints in the first generation 
and the self-adaptive algorithms are all key techniques to improve the search efficiency of Muse. To 
show the search ability of Muse and how it works, I presented some predicted systems including 
metallic, covalent and ionic systems. All the test results show Muse has high efficiency and almost 
100% success rate. 



INTRODUCTION 



The structure of a material provides key information 
for us to conveniently and deeply explore its properties. 
Especially, the first-hand structure information is the be- 
ginning to discover new physics and chemistry of mate- 
rials in extreme environments. Commonly, the structure 
information can be achieved from experiment via the X- 
ray powder diffraction technique. However, when the 
material is in very complex or extreme conditions such 
as high pressure and/or high temperature 1 or its com- 
position is beyond our chemical intuition, 2 experiment 
is often difficult or even impossible to identify its struc- 
ture. On the contrary, the theoretical method is often 
efficient, economic and convenient. Even though Mad- 
dox claimed the "scandal in the physical sciences" of the 
inability to predict the structures of materials from sto- 
ichiometry two decades ago, 3 up to date many groups 
have made great success in structure prediction only from 
stoichiometry based on different stochastic search algo- 
rithms, including the simulated annealing^ the minima 
hopping, 5,6 the basin hopping, 7 the evolutionary algo- 
rithnijii 2 -^ the particle swarm optimization ^ 10 ! 11 and the 
random search ji 2 - Utilizing these methods, many groups 
have found plenty of novel materials under ambient and 
extreme conditions. 

The evolutionary algorithm was inspired by biological 
evolution, such as mutation, reproduction, recombina- 
tion, and selection. The whole population's evolution is 
the process of the survival of the fittest. For the crys- 
tal structure prediction, this algorithm is utilized by in- 
troducing several evolutionary operators. The evolution 
speed is determined by the evolutionary operators in- 
cluding the variation operators and the selection opera- 
tor. The variation operators control the diversity of the 
population, and the selection operator mainly determine 
the direction of evolution. The simulated annealing algo- 



rithm was from annealing in metallurgy, which hastens a 
melted material to be crystal by controlled cooling. Sim- 
ulated annealing improves its efficient through the intro- 
duction of two tricks: the so-called Metropolis selection 
rulei 3 - and lowering the "temperature" . The basin hop- 
ping algorithm is just as its name that means trying to 
escape the local minima by hopping. The basin hopping 
algorithm also uses Metropolis rule^i 3 - 

The optimization algorithms are all stochastic tech- 
niques and have their own inherent advantages and dis- 
advantages. If we can organically combine some of these 
algorithms to let them work collaboratively, we will avoid 
the disadvantages and fully utilize their advantages. In 
Muse, I organically combined the evolutionary, the sim- 
ulated annealing and the basin hopping algorithms to 
collaboratively predict the stable and metastable struc- 
tures of materials. 

The implementation of multi algorithms collaboration 
are based on the manipulation of crystal structures in 
real space, i.e. the operations on crystal structures made 
by the variation operators. So in order to manipulate 
crystal structures conveniently, the structure of crystal is 
described by the shape of the lattice cell and the atoms 
positions therein. The fitness for the survival of the best 
is the free energy of each structure. The structure with 
lowest free energy have the best fitness. Thus, the se- 
lection operators is to select lower free energy structures 
(parents) to generate new structures (offspring). We also 
refer to each structure as an individual and all the struc- 
tures as the population. 

In this paper, I focus on the implementation of the 
multi algorithms collaboration technique and how it 
works. This paper is organized as follows. Section [TT] is 
the algorithms implementation. This includes multi algo- 
rithms collaboration, the descriptions of the evolutionary 
and variation operators, the evaluation and selection of 
an individual and how the search is terminated. I will 
show some results predicted by Muse in Section Hill The 
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predicted systems includes metallic, covalent and ionic 
systems. The discussion and conclusion are provided in 
Section ED 



II. ALGORITHM IMPLEMENTATION 

Muse is developed for easy use in structure predic- 
tion of materials under ambient or extreme conditions, 
e.g. high pressure. It was written in Python and organi- 
cally combined the multi algorithms including the evolu- 
tionary algorithm, the simulated annealing algorithm and 
the basin hopping algorithm to collaboratively search the 
global energy minima of materials with the fixed stoi- 
chiometry. After introduced the competition in all the 
evolutionary and variation operators, the evolution of 
the crystal population and the choice of the operators 
are self-adaptive automatically. That is to say the crys- 
tal population undergoes a self-adaptive evolution pro- 
cess. So, it can effectively find a material's stable and 
metastable structures under certain conditions only pro- 
vided the chemical information of the material. 

Muse depends on external local optimization codes. 
At present, Muse uses VASP,i^£ SIESTA^ Quantum 
ESPRESSO 17 and LAMMPS 18 as the local optimization 
tools. It determines the space group number of each op- 
timized structure immediately after the optimization is 
done. The duplicate structures are then eliminated ac- 
cording to the space group numbers and the enthalpies. 
In detail, the supercells with the same space group num- 
ber are first reduced to primitive cells. Then the dupli- 
cated structures are deleted if the numbers of atoms in 
the primitive are equal and their enthalpy difference are 
less than 1 meV/atom. More importantly, Muse gener- 
ates the random structures of the first generation with 
symmetry constraints, which largely shortens the opti- 
mization time of the first generation and increase the di- 
versity of crystal population. The random structures are 
created according to the randomly chosen space group 
numbers from 2 to 230 and Wyckoff positions must be 
fit to the atom numbers ratio of different kinds. Espe- 
cially for large systems, the constraints on symmetry will 
avoid the unphysical disorder crystal structures, similar 
to glass state. Muse can also pick up (restart) a previous 
interrupted search from wherever the search stopped. 



A. Multi-algorithm collaboration 

A single algorithm has more or less inherent disad- 
vantages. Using the combination of multi algorithms, 
or hybrid, we will avoid their disadvantages and effec- 
tively make use of their advantages. The efficiency of the 
multi-algorithm collaboration search also depends largely 
on the diversity of crystal population, which is up to the 
evolutionary and variation operators. To increase the 
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FIG. 1: The flowchart of multi algorithms collaboration and 
coupling in Muse 



diversity of crystal structures, MUSE uses ten evolution- 
ary and variation operators, including cross over, mu- 
tation, permutation, cross-over-mutation, permutation- 
mutation, slip, twist, random move, ripple and mutation- 
ripple. The multi-algorithm collaboration and coupling 
are shown in Fig. [TJ Different systems have different op- 
timal operators to achieve the possible largest diversity. 
So, from the second generation these ten evolutionary 
operators will compete with each other in their success 
rates starting from the initial value 100% of each opera- 
tor. The operators are randomly used according to their 
historical behaviors from the third generation. That is 
to say, the more positive contributions (produced lower 
enthalpy structures or increased the diversity) an oper- 
ator has done to the whole population previously, the 
more chances it will have to be called to breed offspring. 
This is called the self-adaptation of evolutionary oper- 
ators for different systems, which turns out to largely 
increase the diversity of structures and hasten the con- 
vergence of global search. 

Whether the new optimized offspring in the main 
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evolutionary loop is kept or not is determined by the 
Metropolis algorithm 1 - based on the decreasing temper- 
ature, which actually couples the evolutionary algorithm 
and the simulated annealing algorithm together. Fur- 
thermore, to prevent prematuration, the main evolution- 
ary loop is also coupled with basin hopping algorithm 
which has advantage of global optimization. After the 
main evolutionary loop is converged, Muse will enter into 
the pure basin hopping loop to try to escape the possible 
local minimum. Tests show if we gather the structure 
information from the small systems with small number 
of primitive cells, we will easily convergence the large 
systems who have more primitive cells. So we are en- 
couraged to start from small systems with small number 
of primitive cells. This process is automated in Muse. If 
all the systems are converged, Muse will terminate the 
search and exit. 




B. The first generation 

A unit cell in Muse is described with six parameters, 
a, b, c, a, j3 and 7, and the atoms positions, where a, 
b and c are the lattice vectors and a, f3 and 7 are the 
corresponding angles as normal. The structures of the 
first generation can be randomly created with the sym- 
metry constraints by randomly choosing crystal space 
group number (with Wyckoff positions constraints) from 
2 to 230 and the constraints of the minimum and maxi- 
mum angles (45° and 135°) between lattice vectors. Or, 
if one wants to generate the first generation in the fully 
random manner, Muse will generate the structures with- 
out symmetry constraints but still with the constraints 
on minimum and maximum angels between lattice vec- 
tors. Then its volume is scaled to the coarsely guessed 
one with the constraints of the minimum length of lattice 
vectors and the minimum distance between atoms. If a 
space group number is used once, Muse will try to avoid 
using it for the rest of the random structures in order to 
prevent duplication and increase the diversity. 



C. Evolutionary and variation operators 

The individuals in the second generation are created 
using the evolutionary and variation operators from the 
certain percentage of the number of individuals in the 
first generation. The ten operators fall into two cate- 
gories: the single-parent based operators and the two- 
parent based operators. The mutation, the permuta- 
tion, the random move, the ripple, the slip, the twist, 
the permutation-mutation and the mutation-ripple op- 
erators are all the single-parent based operators. The 
cross over and the cross-over- mutation operators are two- 
parent based operators. 



FIG. 2: The example of mutation, (a) Randomly chosen par- 
ent: P-l (2), (b) Offspring: P2i/m (11). The two structures 
are locally optimized. 



1. Mutation 

The mutation operator variate one selected parent to 
produce a new individual (offspring) by multiplying the 
unit cell row vector with a symmetric Voigt strain ma- 
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where 5{j are the zero-centered Gaussian random vari- 
ables with a specified standard deviation, a is the orig- 
inal lattice vector, and a 7 is the new lattice vector. The 
atomic coordinates are fractional values and are scaled 
to the new vectors. The mutation operator is actually 
applying a shear strain on the cell to cause it to undergo 
phase transition. So the mutation operator is the very 
effective operation to diversify the crystal population. 
However, different systems have different optimal stan- 
dard deviations 5{j, so Muse uses trial and error scheme 
to obtain good diversity from the initial standard devi- 
ation value. After applying the strains the crystal's vol- 
ume is then scaled to the volume of the lowest enthalpy 
structure in the previous generation. For the example of 
mutation, please see Fig. [2j 



2. Permutation 

The permutation operator changes the positions of two 
kinds of atoms a random number of times. 1,2,19 This ob- 
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viously increases the diversity of the crystal population 
after operation. The permutation operator can only be 
called when the atom types is more than one. 

3. Permutation-mutation 

The permutation-mutation operator is the hybrid of 
the permutation and mutation operators. The unit cell 
is mutated after exchanging two atoms of two different 
kinds random number of times. It is the single-parent 
operator, but it obviously increases the diversity of the 
crystal population. 

4- Random move 



the coordinates of each atom along a randomly chosen 
axis by certain amount. For example, if x-axis is chosen, 
the x components of all atoms will be shifted to the new 
values in fraction: where Ax depends on 

the atoms non-displaced (i.e. y and z) coordinates via* 2 - 

Ax — pcos(27Tjj J y + 0y)cos(27rr]z + 6 Z ), (2) 

where p is the maximum possible displacement in x direc- 
tion, and fx and 77 are integers, and y and 6 Z are random 
numbers between [0, 2ir). 

The ripple operator is also a singe-parent operator, but 
it can substantially increase the crystal diversity. Some 
solids manifest the ripple motif at ambient or high pres- 
sure, such as Cs-III, 20 Rb-III 21 and Gn-IlM Fig. [3] is an 
example of ripple. 



The random move operator only randomly moves the 
atoms positions by random amount, keeping the lattice 
vectors unchanged. It displaces the three fractional coor- 
dinates of all atoms by random amount between [-D max , 
D max ], where D max is the maximum percentage of co- 
ordinates that atoms displaced. The use of this opera- 
tor combined with the Metropolis algorithm realized the 
basin hopping loop^^ The periodic boundary conditions 
are applied after randomly moved atoms. This loop is 
coupled with the main evolutionary loop. 



5. Slip 

To fully variate the crystal structure, Muse introduced 
a new operator slip. The slip operator let a group of 
atoms slip by a random distance along a random direction 
that is parallel to a specific plane. This operator comes 
from the real crystal phase transition process in which 
some atom layers or atoms slip along special direction. 
So it will considerably increase the diversity of crystal 
structure. 



6. Twist 
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FIG. 3: The example of ripple, (a) Randomly chosen parent: 
PI (1), (b) Offspring: Cmcm (63). The two structures are 
locally optimized. 



The twist operator is also a new operator first con- 
structed in Muse. It twists the supercell by rotating 
atoms along the randomly chosen axis. After this opera- 
tion atoms will be rotated by #0 + A# around the chosen 
axis, where A# is the interval of rotated angles between 
different atoms. Twist often causes phase transition in 
real crystal under extreme conditions. So it also a high 
efficiency operator. 



8. Mutation-ripple 

This operator is the hybrid of the mutation operator 
and the ripple operator. The unit cell is acted by the 
ripple operator after it is mutated. 



9. Cross over 

7. Ripple 

The cross over operator is the two-parent based opera- 
The ripple operator, first proposed by Lonie and tor. It cuts and splices the randomly picked two parents 
Zurekj 2 - is the periodic displacement operator which shifts to make a new offspring 1 9 1 23 ~ ^£ It mainly consists of 
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two steps: cut and splice. In the first step, the cut plane 
is parallel to the randomly not picked two vectors of one 
randomly picked parent and the cut position is in the 0.5- 
centered Gaussian random number between [0.25, 0.75]. 
If the atoms coordinates of the picked axis is below the 
Gaussian random number, they are taken as the part of 
the new individual. The other part of the new individual 
is from the other parent in the same way but the atoms 
coordinates are larger than the same Gaussian random 
number. In the second step, the picked atoms from the 
two parents are placed together to make the atoms coor- 
dinates of new individuals. The new lattice vectors are 
the vectors summation of the two parents. To increase 
the success rate and effectiveness, the atoms are centered 
before cross over operation. Fig.|4]is one example of cross 
over. 



determined by 
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FIG. 4: The example of cross over, (a) Randomly chosen 
parent 1: Ccmm (63), (b) Randomly chosen parent 2: P6s/m 
(176), (c) Offspring: Pm (6). The three structures are all 
locally optimized. 



Hi 



(3) 



where Hi is the enthalpy of this individual, and H max 
and i^min are the enthalpies of the worst and the best 
individuals in the last generation, respectively. 

After an offspring is generated and locally optimized, 
whether it is kept or not is up to the Metropolis rule^: 

• An individual is kept if its enthalpy is lower than 
that /those of its parent (s). 

• But if its enthalpy is greater than that/those of its 
parent (s), it is kept only if 



GXp [(-^parents -^individual ) 

is greater than a randomly picked number from [0, 1]. k 
is Boltzmann constant. The initial temperature Tq for 
annealing is specified in the input file. The temperature 
in generation n is reduced to T n = 0.9T n _i. 



E. Termination 

The termination of the main evolutionary loop coupled 
with the simulated annealing and the basin hopping in 
Muse is controlled by the terminator operator. The main 
loop is terminated: 

• if the number of continuous generations with the best 
enthalpy difference less than 1 meV/atom and the same 
symmetry structures reaches the number specified in the 
input file; 

• if the specified maximum number of generations has 
been done; 

• if the diversity of the crystal population is too low. 
While, the pure basin hopping loop is terminated when 

the loop number reaches the number of continuous gener- 
ations with the best enthalpy and the same space group 
number. 



10. Cross-over-mutation 

The cross-over-mutation operator is the hybrid of the 
cross over operator and the mutation operator. 



D. Evaluation and selection 

As other stochastic structure search algo- 
rithmsji»2&i2r— Muse also uses ab initio Gibbs 
free energy as the fitness function, which evaluates if an 
individual is suited to be parent to breed offspring. At 
K, Gibbs free energy equals to enthalpy. The enthalpy 
(H = E + PV) of an individual is determined from the 
local optimization code. The absolute probability pi of 
an individual being selected for breeding offspring is 



F. K-point adaptation 

The K-point adaptation is similar to Ref. yj. The K- 
grid in reciprocal space of a lattice vector i is calculated 
via: 



Oj\ • rC reso l 

where a\ is the length of the lattice vector z, and k reso \ is 
the reciprocal-space resolution set in the input file. k[ is 
then rounded to an integer. 

III. RESULTS 

The multi-algorithm collaboration in Muse have been 
sufficiently tested for many cases including metallic, co- 
valent and ionic systems. Table [I] shows the searched 
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TABLE I: The stable structures of different systems predicted 
by Muse under different conditions. N a tom is the number of 
atoms in the system. P is the target pressure in GPa. Size is 
the population size in each generation. Nq shows how many 
generations Muse used to find the stable structure. 



System A^tom 


P (GPa) 


Population size 


Symmetry Nq 


Pd 


10 





20 


Fm3m 1 


Ta 


10 





20 


Im3m 1 


Au 


8 





20 


Fm3m 1 


Nb 


8 


100 


20 


Im3m 1 


Y 


8 


140 


30 


Fddd 2 


Sr 


20 


100 


30 


C2/m 3 


C 


12 


50 


30 


Fd3m 4 


SiC 


20 


20 


30 


F43m 2 


LiBC 


12 





30 


P63/mmc 6 


NaCl 


20 





30 


Fm3m 1 


GaN 


20 





30 


P6 3 mc 2 


NbN 


8 


100 


30 


P63/mmc 1 


NbS 2 


18 


50 


30 


I A 1 rnrnrn 3 


MgSi0 3 


20 


130 


30 


Cmcm 5 



stable structures of different systems under different con- 
ditions. The random structures in the first generation 
were constructed with symmetry constraints. The ab ini- 
tio optimizations and the free energy calculations for ev- 
ery structure generated by Muse were performed with 
VASP.i^ The GGA parametrized by PBE 29 was ap- 
plied and the electron-ion interactions was described by 
the PAW scheme . 3Q i 31 To achieve good convergences the 
kinetic energy cutoff and the Appoint grids spacing were 
chosen to be 1.3 times the default values and 0.03 A , 
respectively. As one can see, the systems with not more 
than two kinds of atoms are very easy for Muse to find 
the stable structures; it generally finds them in not more 
than three generations. However, the systems with more 
than two kinds of atoms need at least five generations. 
On the other hand, in accord with our intuition larger 
systems take more generations to successfully find the 
global lowest enthalpy structures, because the configura- 
tion space increases exponentially with the system size. 
All the searched stable structures are in agreement with 
the known structures. 



A. Metallic systems 

Tests on metallic systems have been done for Pd, Ta, 
Au, Y, Nb and Sr. Except for Y and Sr, the known 
stable structures of other metals were found in the first 
generation. Y is a metal at ambient conditions. Un- 
der compression, it will exhibit superconducting prop- 
erties according to recent report^ At 140 GPa, Muse 
also found the superconducting Fddd phase (see Table U 
for details) recently reported by Chen et. al. 32 Apart 
from the Fddd phase, Muse also found other energeti- 
cally competitive structures, such as the structures with 




(b) 



FIG. 5: The predicted structures for 12 C atoms at 50 GPa. 
(a) graphite, (b) diamond. 



C2/c (No. 15) and P2\jc (No. 14) which have not been 
reported before. 

B. Covalent systems 

The covalent systems including C, SiC and LiBC were 
also fully tested. After performed search for 12 C atoms 
at 50 GPa, Muse successfully reproduced the diamond- 
C and graphite-C in the same search (see Fig. [5j). The 
diamond phase is the most stable one. Other metastable 
structures with symmetries Cmmm (No. 65), Fmmm 
(No. 69), R3m (No. 166), C2/m (No. 12), P62ra 
(No. 189) and so on, are also found. SiC with two atom 
types only takes three generations for Muse to find the 
stable structure with space group Fm3m. Many ener- 
getically competitive structures are also found, including 
structures with P6smc, Ccra2i, C2/ra, Cm and P3ml 
symmetries. The enthalpy of the structure with P6smc 
(No. 186) symmetry is only 3.1 meV/atom higher than 
that of Fm3m at 10 GPa. LiBC is an intermetallic com- 
pound which also shows strong in-plane covalent bond- 
ing. Muse only used 6 generations to find the known 
P6s/mmc structure for the 12 atoms system. 

C. Ionic systems 

NaCl is a typical well-known ionic system. With the 
help of symmetry constraints on the first generation, 
Muse found its stable Fm3m in the first generation using 
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012345678 
Generation Number 

FIG. 6: The self-adaptive evolution of SiC structures at 10 
GPa. The cell contains 10 atoms. The structure with F43m 
symmetry was found in the third generation. The inset shows 
all the structures searched. 




Cmcm Cmcm 



2 4 6 8 10 

Generation Number 

FIG. 7: The self-adaptive evolution of MgSiOs structures at 
130 GPa. The cell contains 20 atoms. The well-known post- 
perovskite structure (Cmcm) of MgSiOs was found in the fifth 
generation. The inset shows all the structures searched. 



20 atoms. Fig. [71 shows the evolution of MgSiOs struc- 
tures at 130 GPa. There are 20 atoms in the unit cell. 



With so many atoms in the cell, Muse find the well- 
known post-perovskite structure only in five generations. 
At the same time, Muse produced many metastable sta- 
ble structures with relatively lower enthalpies. Smaller 
MgSiOs systems with 10 atoms were also tested at 130 
GPa, and Muse only used two generations to find Cmcm 
structure. 

IV. DISCUSSION AND CONCLUSION 

The implementation of the multi-algorithm collabora- 
tion for crystal structure prediction are detailed. After 
introduced two new variation operators, slip and twist, 
the search ability of Muse are enhanced. In order to fur- 
ther increase the search efficiency of Muse, the number 
of variation operators are also increased to ten. Further- 
more, after adopting the competition scheme in the ten 
variation operators, different systems will find different 
optimal variation operators in evolution and then fas- 
ten the convergence speed. These techniques are all key 
to improve the efficiency of crystal structure prediction. 
Tests on different systems show the very high efficiency 
of multi-algorithm collaboration. 

Using 12 C atoms, Muse successfully reproduced the 
known diamond and graphite structures, apart from some 
energetically competitive metastable structures. As for 
the complex MgSiOs systems with 20 atoms, Muse 
only used five generations to find its well-known post- 
perovskite phase at 130 GPa. For the metal systems, 
tests on Ta, Pd, Au, Y, Nb, Sr, etc, also show very high 
efficiency. For metal Y, Muse found the Fddd phase at 
100 GPa, which was recently reported. In all tests, Muse 
achieved the success rate of almost 100%. So, the multi- 
algorithm collaborative crystal structure prediction has 
high search efficiency and success rate. 
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